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We performed two-dimensional simulated tempering (ST) simulations of the two-dimensional Ising 
model with different lattice sizes in order to investigate the two-dimensional ST's applicability to 
dealing with phase transitions and to study the crossover of critical scaling behavior. The external 
field, as well as the temperature, was treated as a dynamical variable updated during the simulations. 
Thus, this simulation can be referred to as "Simulated Tempering and Magnetizing (STM)." We 
also performed the "Simulated Magnetizing" (SM) simulations, in which the external field was 
considered as a dynamical variable and temperature was not. As has been discussed by previous 
studies, the ST method is not always compatible with first-order phase transitions. This is also 
true in the magnetizing process. Flipping of the entire magnetization did not occur in the SM 
simulations under T c in large lattice-size simulations. However, the phase changed through the 
high temperature region in the STM simulations. Thus, the dimensional extension let us eliminate 
the difficulty of the first-order phase transitions and study wide area of the phase space. We then 
discuss how frequently parameter-updating attempts should be made for optimal convergence. The 
results favor frequent attempts. We finally study the crossover behavior of the phase transitions 
with respect to the temperature and external field. The crossover behavior was clearly observed in 
the simulations in agreement with the theoretical implications. 

PACS numbers: 64.60.De,75.30.Kz,75.10.Hk,05.10.Ln 

Keywords: multi-dimensional generalized-ensemble algorithm, two-dimensional simulated tempering (ST), 
simulated tempering and magnetizing (STM), Monte Carlo (MC) simulation, Weighted Histogram Analysis 
Method (WHAM), Multistate Bennett Acceptance Ratio estimator (MBAR), Ising model, phase transition, 
critical phenomenon, crossover 



I. INTRODUCTION 

In the computational statistical physics field, Monte 
Carlo (MC) and molecular dynamics (MD) simulations 
have been commonly used. However, the quasi-ergodicity 
problem, where simulations tend to get trapped in states 
of energy local-minimum, has often posed a great diffi- 
culty. In order to overcome this difficulty, generalized- 
ensemble algorithms have been developed and applied to 
many systems including spin systems and biomolecular 
systems (for reviews, see, e.g., Refs. [1-3]). 

Commonly used examples of generalized-ensemble al- 
gorithms are multicanonical (MUCA) [4, 5], simulated 
tempering (ST) method [6, 7], and replica-exchange 
method (REM) [8, 9] (it is also referred to as parallel 
tempering). Closely related to MUCA are the Wang- 
Landau method [10, 11] and metadynamics [12]. Also 
closely related to REM is the method in Ref. [13]. 

In the ST method, temperature is regarded as a dy- 
namical variable, which is updated by the Metropolis cri- 
teria during the simulation, and consequently a random 
walk is realized in the temperature space. This random 
walk, in turn, causes a random walk of the energy, which 
enables the system in question to overcome free-energy 
barriers. However, it is well-known that the ST method 
is not very compatible with first-order phase transitions 



(for a review, see, e.g., Ref. [14]). When there is a first- 
order phase transition, the random walk of temperature 
across the phase-transition point hardly occurs. We re- 
mark that there is a recent attempt to deal with this 
difficulty by an extension of ST [15]. 

Recently, the multi-dimensional generalizations of the 
generalized-ensemble algorithms, including MUCA, ST, 
and REM, were discussed and general formalisms were 
given [16-18]. In these methods, the energy of the system 
is generalized by adding other energy term(s) with some 
coupling constants. In the multi-dimensional ST method, 
not only the temperature but also the coupling constants 
are considered as dynamical variables. 

In this work, we study a special case of the above gen- 
eral multi-dimensional ST methods. Namely, the addi- 
tional term is —hM where h and M are the external field 
and the magnetization, respectively. The external field 
h corresponds to the couping constant which is updated 
during MC simulations. Therefore, not only temperature 
but also external field becomes a dynamical variable and 
is expected to realize a random walk during the simu- 
lations. Thus, this simulation can be referred to as the 
"Simulated Tempering and Magnetizing" (STM). In or- 
der to test the effectiveness of the present method, we 
applied it to the two-dimensional Ising model. 

The Ising model has two kinds of phase transitions. 
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One occurs along the change of temperature when the ex- 
ternal field is zero. The other occurs along the change of 
external field when the temperature is under the critical 
temperature (T c ). The former is classified into a second- 
order phase transition. The latter is categorized as a 
first-order phase transition unless the temperature is ex- 
actly equal to T c . When T = T c , the transitions are clas- 
sified into a second-order phase transition. This system 
allows us to confirm applicability of the two-dimensional 
ST to the first-order phase transitions along the external 
field changes. 

We also investigated the crossover phenomena in 
the phase transitions, in which critical exponents are 
changed. We study the behavior of magnetization per 
spin m, which follows m ~ \T — T c \@ and m ~ j/il 1 /" 5 near 
the critical point, where /? and S are critical exponents 
[19]. Our simulation method, with a combination of his- 
togram reweighting techniques, enables us to calculate 
physical values such as the energy and magnetization at 
various values of T and h from a single production run. 

This article is organized as follows. In §2 we present 
the STM method. In §3 we present the results. In §4 we 
conclude this article. 



II. MATERIALS AND METHODS 
A. System 

We study the two-dimensional Ising model in external 
field. The total energy is given by 

H = E — hM , (1) 
E = - , (2) 

(id) 

N 

M = J2*i, (3) 

i=l 

where i, N, cri, and h are the index of spin, total number 
of spins, spin at the i-th site, and external field, respec- 
tively. The spin cr, takes on the values ±1. The sum in 
Eq. (2) goes over the nearest-neighbor pairs. The spins 
are arranged on the square L x L lattice. We imposed 
the periodic boundary conditions. Data were obtained 
for lattice sizes from 2 x 2 to 160 x 160. 



B. Simulation methods 

Whereas the conventional ST method considers tem- 
perature as a dynamical variable, the STM method con- 
siders not only temperature but also external field as a 
dynamical variable. Here, before explaining the STM 
method, we shortly review the conventional ST method 
[6, 7]. 

In the conventional ST method, temperature is a dy- 
namical variable which takes on one of Nt values (here, 



temperature is discretized into Nt values). In other 
words, denoting X and x as a sampling space and its 
microscopic state, respectively, the Boltzmann factor 

e -E(x)/T+a(T) ^ 

is regarded as a joint probability for the state (x, T) 
{Ti,T 2 , . . . , Tjv t }, the product space of X and 
{T\,T2, ■ ■ ■ ,Tat t }). Here, a(T) (or a(Tj)) is a parameter 
for obtaining uniform distributions of temperature val- 
ues. Here and hereafter, we set Boltzmann's constant to 
unity. Now that the temperature is a dynamical vari- 
able, the simulated system is allowed to realize a random 
walk in the temperature space. This random walk, in 
turn, causes a random walk of energy. Consequently the 
simulated system has more chance to overcome energy 
barriers. 

Even though temperature changes during ST simula- 
tions, any thermodynamic quantity at temperature T,, 
(A) T ., can be reconstructed with the conditional expec- 
tation of a physical quantity A given at T,, or (A\Ti). 
Note that 



(A\Ti) ST — 




where Sij is the Kronecker delta. Namely, we have 

(8) 

" 3 = 1 

where N^ and A T . stand for the total number of samples 
and j-th sample at Ti. 

To find the candidate for a(Tj), let us look into the 
probability of visiting T^. By summing over the delta 
function, the probability of occupying Tj is given by 
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where / is the dimcnsionless (Hclmholtz) free energy and 

= J dxe- E ^' T . (12) 

Substituting /(Ti) into a{Ti) gives constant probability 
regardless of Ti. Thus, the dimcnsionless free energy 
f(Ti) is a good choice for a(Ti) in order to obtain uni- 
form temperature distribution and to realize a random 
walk in the temperature space. Although the free en- 
ergy is not known a priori, unless the system is exactly 
solvable, the free energy calculation methods (the details 
will be provided below) enable us to get its good estimate 
from preliminary simulation runs. 



In the two-dimensional ST algorithm, on the other 
hand, we consider that another parameter is also a dy- 
namical variable [16-18]. Especially in the STM method, 
the external field h is a second dynamical variable. In 
other words, we consider 

e -(E-hM)/T+a(T,h) / 13 N 

as a joint probability for (x,T,h) (e X ® 
{Ti,T 2 ,...,T Nt } ® {hi,h2,...h Nh }), where a(T,h) 
is a parameter. 

To find the candidate for a(Tj, hj), we again look into 
the probability of staying at each set of parameter values. 
It is given by 



N T N h 



P{Ti,hj 



k=l 1=1 



E{x)-h l M{x) 



+a(T k ,h,) 



N T N h 



dxe 

k=l 1=1 J 

N T N h 



. EM -h,M(*) +a(Tfc;M 



oc e 



fe=i i=i 



f{T k ,h l )+a(.T k ,h i ) 



(14) 



(15) 



(16) 



where 

e -m, hj ) = J dxe -{E- hj M)/T K (17) 

The dimcnsionless free energy f{Ti,hj) is again a good 
choice for a(Ti,hj) in order to acquire uniform distribu- 
tion of T and h. These values can be estimated from 
preliminary simulation runs and reweighting techniques. 

As in conventional ST method, any thermal aver- 
age {A) Tuh . at given T!; (e {T x , T 2 , . . . , T Nt }) and h 3 
(g {hi, h%, . . . , hpf h }) can be obtained by calculating 



the conditional expectation: (A) T . h , = (A\Ti, hj) ST . 
Namely, we have 

(a)t /j = j^— A khj . ( 18 ) 

where Nx^hj is the total number of samples at Ti and hj, 
and At^. h . stands for the k-th sample at Tj and hj . 

The method of updating T or h is similar to that of 
updating spins because T and h are considered as dy- 
namical variables. The Metropolis criterion for updating 
T or h is given by the following transition probability: 



^hj^T^hj^rn^^M 



(19) 



Once an initial state is given, the STM simulations can 
be performed by repeating the following two steps. 1. We 
perform a conventional canonical simulation at Ti and hj 



for certain MC sweeps. 2. We update the temperature or 
external field by Eq. (20) with a(T, h) = f(T, h). 

In our implementation every certain MC sweeps ei- 
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ther T or ft, was updated (the choice between T and h 
was made at random) by Eq. (20) to a neighboring value 
(the choice of two neighbors was also made at random). 
Here, one MC sweep stands for Lx L single spin updates. 
The number of MC sweeps performed between parame- 
ter updates is here referred to as the parameter-updating 
period. 

Whereas updating the parameter to a neighboring 
value with the Metropolis algorithm should be consid- 
ered the easiest to implement, we remark that, as spins 
can be updated by a number of methods such as the heat 
bath method, other schemes of updating the parameters 
can be employed [20]. There also exists a temperature 
updating scheme for ST by Langevin algorithm [21]. 

Table I summarizes the conditions of the present sim- 
ulations. For L = 80, instead of a single 4000000000 MC 
sweep production run, four 1000000000 MC sweep runs 
were performed. This was just to make one trajectory 
shorter and easier to deal with numerically. Similarly, 
two production runs (instead of a single run) were made 
for L = 30 and 160. 

As for spin-updates, we employed the single spin up- 



n(E, M) 



where riT u h (E,M) is the histogram of E and M at Tj 
and hj, and Nx^hj is the total number of samples ob- 
tained at Ti and hj . By solving these two equations self- 
consistcntly by iterations, we can obtain n(E, M) and 
f(Ti,hj). The obtained n(E,M) allows one to calculate 
any thermal average at arbitrary temperature and exter- 
nal field values. Note that /(Tj, hj) are determined up to 



where N, NT k .h n E n , and M n is the total number of 
data, the number of samples associated with and hi, 
energy of the n-th data, and magnetization of the n- 
th data, respectively. This equation should be solved 
self-consistently for f(T t ,hj). Note that, as in WHAM, 



date algorithm; we updated spins one by one with the 
Metropolis criteria. As for quasi-random-number gener- 
ator, we used the Mersenne Twister [22]. 



C. Free energy calculations 

The simulated tempering parameters, or free energy 
in Eqs. (13) and (17) can be simply obtained by the 
reweighting techniques applied to the results of prelim- 
inary simulation runs [16-18, 23]. We employed two 
reweighting methods for this free energy calculation. One 
method is the multiple- histogram reweighting method, or 
Weighted Histogram Analysis Method (WHAM) [24, 25] 
and the other is Multistate Bennett Acceptance Ratio 
estimator (MBAR) [26], which is based on WHAM. 

The equations of WHAM algorithm applied to the sys- 
tem is as follows. For details, the reader is referred to 
Refs. [17, 25]. The density of states (DOS) n(E, M) and 
free energy values /(Tj, hj) can be obtained by 



(21) 



(22) 



I 

a constant, which set the zero point of free energy. Ac- 
cordingly, n{E, M) are determined up to a normalization 
constant. 

The MBAR is based on the following equations. 
Namely, by combing Eqs. (21) and (22), the free energy 
can be written as 



(23) 



I 

f(Ti, hj) are determined up to a constant. 

We repeat the preliminary STM simulations and free 
energy calculations until we finally obtain sufficiently ac- 
curate free energy values which let the system perform a 
random walk in the temperature and external field space 



£ n Ti , hj {E,M) 

= 

" J2 N Ti,hj exp(/(Ti, hj) -(E- hjM)/Ti) ' 

Ti,hj 

= - log J2 n ( E > M ) cxp(-(£ - hjM)/Ti) , 

E.M 



t/m u \ i V^V^V^ exp(-(£„ - hjM n )/Tj) 

f(T h hj) = - io g ^ ^ X, -wir h 

i-i 3 -i n-i ^2 N Tk , hl exp(/(T fc , hi) - {E n - hiM n )/T k ) 



fe=i i=i 
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TABLE I. Conditions of the two-dimensional ST simulations. 



Lattice size L 2, 4, 8, 10, 20 

Number of production runs 1 

Total MC sweeps per run 42000000 

Parameter-updating period 50 

T 1 -T Nt 1.0-5.0 

hi-tiN h -1.5-1.5 

N T 20 

N h 21 

iVdata" 10 



a The data were stored every -/V<j ata MC sweeps. 

during the STM simulation. We then perform a single, 
final production run. 

Note that these two reweighting methods enable us to 

I 

<A)t,h 



(Ca) 



For details, the reader is referred to Refs. [26, 27]. 

We also used another method of calculating free en- 
ergy. By substituting a(T, h) in Eq. (16) by the estimates 
for free energy f(T, h), we obtain 

P(T,h)oce-^ T - h ^ T ' h l (27) 

From this we can write 

f(T, h) = f(T, h) - log P(T, h) + const. (28) 

Here, P{T, h) can be obtained as the number of samples 
at each set of parameter values in a preliminary STM 
simulation. Thus, this equation enables one to refine the 
free energy much more easily than the reweighting meth- 
ods, because the method does not require any iterations. 
This method does not work well, however, when P(T,, hj) 
is too small (or /(Tj, hj) is too far away from true values) 
to obtain samples at (Tj, hj), while the reweighting tech- 
niques are still able to work. In the present work, we first 
used the reweighting methods to obtain rough estimates 
of the free energy for the entire parameter space. We 
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obtain not only dimensionless free energy values but also 
physical values at any temperature and at any external 
field. It is given by 



(24) 
(25) 



(26) 



I 

then used the combination of the reweighting methods 
and Eq. (28) for further refinements of the free energy. 



Note that the WHAM gives another piece of informa- 
tion, namely DOS, which MBAR cannot directly calcu- 
late. However, the WHAM requires to make histograms 
before iterations and two kinds of calculations in an it- 
eration step. As the system size grows, the number of 
possible states increases. Thus, the calculation of DOS 
can be quite time-consuming. On the other hand, MBAR 
can be used without making histograms and one MBAR 
iteration step needs one equation. The length of one 
iteration, which is approximately proportional to the 
number of samples and parameter values, increases and 
can be time-consuming, as the system size is enlarged. 
However, we have an impression that the MBAR is less 
time-consuming and more easily implemented than the 
WHAM. The parallelization of MBAR is slightly eas- 
ier than that of WHAM and we actually did it with 
OPENMP. 



E W na A{x n ) , 



J exp{-(E n - hM n )/T) 

(c a ) N ^ N » 

E E N ^M exp(/(T fe! hi) - {E n - h[M n )/T k ) 

k=l 1=1 

eM-(E n -hM n )/T) 

" =1 E E N n,h, cx P (/(T fe , hi) - {En - hiM n )/T k ) 

k=l 1=1 
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D. Temperature and external field distributions 

As is mentioned in the previous subsections, we have to 
give the set of temperature and external field values be- 
fore ST or STM simulations. Actually the determination 
involves trial and error. However, still the reweighting 
methods help one to do this to a certain extent. 

Firstly the maximum and minimum values of temper- 
ature and external field were chosen so that the area of 
temperature and external field were wide enough to in- 
vestigate the critical behaviors. This should be done sep- 
arately for each system and what is to be investigated. 

The distribution of temperature was chosen to be pro- 
portional to an exponential to the index number in small 
lattice sizes, as is common in simulated tempering and 
replica-exchange methods. However, in large lattice size 
systems, we assigned more number of values around T c 
by hand. More dense distribution was required where 
the heat capacity is large or the phase transition occurs. 
The distribution of external field is similarly assigned. In 
small lattice size it was proportional to the index of exter- 
nal field. However, in the larger lattice size, we assigned 
more points around h = 0, in which the phase transition 
occurs. We assigned them in such a manner that the 
acceptance rate of ST parameter updates arc preferably 
between 10% and 50%. This fuzzy criterion is partly due 
to the two-dimensional distributions. A temperature dis- 
tribution at a certain external field does not always give 
the same acceptance rates under another external field. 

When the distributions of Tj and hj turned out to be 
improper, we reassigned the distributions. In this case, 
we already had the samples and free energy estimates 
at a previous distribution, with which the reweighting 
method lets one to estimate the free energy at the newly 
distributed values. Consequently, we did not have to 
start over the free energy calculations from the begin- 
ning. We actually repeated this parameter redistribution 
procedures several times, especially in large lattice size 
simulations. 



III. RESULTS AND DISCUSSION 

A. "Simulated Tempering and Magnetizing" 
(STM) simulations 

Firstly we shall show that the two-dimensional ST sim- 
ulations were carried out properly. Figures 1 and 2 show 
temperature and external field, respectively, as functions 
of MC sweep. Both were obtained from the simulations 
in which the linear lattice size was 80. The temperature 
and external field indeed realized random walks. 

Figures 3 and 4 show energy and magnetization per 
spin, respectively, as functions of MC sweep. They also 
realized random walks. Note that there are expected cor- 
relations between the temperature and energy (see Figs. 
1 and 3) and between the external field and magnetiza- 
tion (see Figs. 2 and 4). The same behavior was observed 




0.5 1 1.5 

MC sweep (x 10 6 ) 



FIG. 1. (Color online) The history of temperature T. The 
linear lattice size L was 80. 




FIG. 2. (Color online) The history of external field h. The 
linear lattice size L was 80. 



in other lattice size simulations (data not shown). 

Figure 5 shows the dimensionlcss free energy per spin 
as a function of temperature and external field, which 
was obtained by applying MBAR to the results of the 
production runs. Note that the partial differential of this 
free energy by h gives The shape at h = suggests 
a jump of m below T c , indicating existence of the first- 
order phase transitions. 

Figure 6 shows the distribution of magnetization as a 
function of temperature. Below T c the distribution is 
separated into two parts. As temperature increases, the 
distribution becomes broader. Near T c the distribution 
is the broadest and two peaks merge. It then becomes 
narrower. Note that this figure was obtained by only four 
production runs (see Table I), and can be obtained even 
by only one production run, though the error is expected 
to become larger. Figures 7(a), 7(b), and 7(c) show the 
distribution of magnetization as a function of external 
field above, around, and below T c , respectively. Above 
T c , the change is smooth and continuous (see Fig. 7(a)). 
Around T c , the distribution becomes very wide around 
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FIG. 5. (Color online) The free energy per spin / /L 2 and its 
contour curves as a function of T and h.The linear lattice size 
L = 80. 



FIG. 7. The distribution of m as a function of h when (a) 
T = 3.21, (b) T = 2.316 and (c) T = 1.967. The linear lattice 
size L was 80. 



h = (see Fig. 7(b)). This is one of the properties of the 
second-order phase transitions. Below T Cl the distribu- 
tion jumps from one side to the other side at h = (see 
Fig. 7(c)). This abrupt jump of distribution is one of the 
properties of the first-order phase transitions. 

We also calculated the Binder cumulant [28] defined 



by 

W)5 5( 3 -^)- (29) 

Figure 8 shows the Binder cumulant as a function of tem- 
perature. As is well-known, the graphs cross at one point 
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at T c . The error bars were obtained by the jackknife 
method [29, 30]. 
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FIG. 8. (Color online) Binder cumulant U vs temperature. 



Figure 9 shows the Binder cumulant MS (X function of 
temperature under different external fields. The graphs 
do not cross at one point in the presence of finite external 
field. The amount of errors is expected to be on the same 
level of Fig. 8 and the error bars are suppressed here to 
aid the eye. 
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FIG. 9. (Color online) Binder cumulant U vs temperature 
under different external fields, (a) h = 0. (b) h = 0.01. (c) 
h = 0.05. (d) h = 0.1. The dotted green curve, solid red 
curve, and dashed blue curve stand for the results for L = 20, 
80, and 160, respectively. 



tional ST shows artifact jumps at a high temperature and 
a certain external field. This must have been caused by 
a failure of sampling some parts of states. On the other 
hand, the results by the STM simulations are smooth 
(see Fig. 11). Figure 12 shows the density of states ob- 
tained by conventional ST and STM simulations. This 
obviously illustrates that the area in which the energy 
is relatively high with somewhat strong magnetizations 
were not sampled by the conventional ST method. These 
results imply that the dimensional extension in the STM 
enlarged the sampled space. 



C. Simulated Magnetizing (SM) 

We study the compatibility of ST with the first-order 
phase transition along external field changes, by perform- 
ing "Simulated Magnetizing" (SM) simulations, in which 
the temperature is fixed and the external field is updated 
by the Metropolis criteria. Figure 13 shows the external 
field as a function of MC sweep in the SM simulations 
below T c . We performed SM simulations in a number 
of lattice sizes from 2x2 to 20x20. These graphs illus- 
trate the fact that as the system size becomes larger, the 
difficulty in simulations grows. In fact it finally became 
impossible to observe the events in which the magneti- 
zation goes to the other side across the zero point (see 
Fig. 14(a)), while it was still possible above T c (see Fig. 
14(b)). These results imply that the full range random 
walk happens above T c but not below T c . Therefore, this 
result suggests that the random walk of temperature is 
crucial for the full range random walk of external field. 
The full range random walk of the external field happens 
in the STM simulation when the temperature was high 
above T c . Note that the Ising model is equivalent to the 
lattice gas model [31]. Hence, what happened in STM 
simulations can be understood as that even though the 
phase transitions between gas and liquid do not directly 
occur, they do occur through the "super critical water 
region." 

To explore this phenomenon more clearly, we present 
Video 1, which shows how the temperature and external 
field changed during the STM simulations. The picture 
shown is the first frame of the video. The red line is 
drawn where the first-order phase transitions occur. The 
movement of parameter can touch the red line but rarely 
go across the red line below T c . On the contrary, above 
T c , the external field can change smoothly. 



B. Comparison of ST with STM 

We compared the results of the STM method with 
those of the conventional ST method. Figures 10 and 11 
show the magnetization as a function of temperature and 
external field, which was calculated using MBAR with 
the data obtained by the conventional ST and STM sim- 
ulation, respectively. Figure 10 obtained by the conven- 



D. How often temperature or external field should 
be updated? 

A common question about this kind of simulation is 
how frequently the parameter-updating attempts should 
be made. We want to emphasize that as long as the 
detailed balance condition is satisfied the simulations 
should be correctly carried out. 
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.Relighted (MBAR) 




FIG. 10. (Color online) Reweighted data (red) and original 
data (green) obtained by the conventional ST. The linear lat- 
tice size L was 80. 



2dim ST 




FIG. 11. (Color online) Reweighted data (red) and original 
data (green) obtained by the STM simulations. The linear 
lattice size L was 80. 



(a) (b> 




FIG. 12. (Color online) Calculated DOS obtained by WHAM 
with ST (a) and STM (b) data, respectively. The linear lattice 
size L was 80. 
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FIG. 13. (Color online) External field vs MC sweep in SM 
simulations under T c (T — 1.97). The linear lattice size L is 
(a) 2, (b) 4, (c) 8, and (d) 10, respectively. 
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FIG. 14. (Color online) External field and MC sweep in the 
SM simulation (a) under T c (T = 1.97) and (b) above T c 
(T = 3.88). The linear lattice size L is 20. 



We compared STM simulations performed with differ- 
ent parameter-updating frequencies. Figure 15 shows the 
results of the heat capacity as a function of temperature 
at h = 0, which were obtained by the STM method with 
different conditions. The conditions are one parameter- 
updating attempt every one, two, twenty, and a hundred 
MC sweeps. They show good agreement with each other. 
The error bars were obtained by the jackknife method 
[29, 30]. Note that the error bars tend to be larger as the 
parameter-updating frequency becomes less. 

Figure 16 shows the magnetization as a function of 
temperature at h = 0. Data were obtained with several 
parameter-updating frequencies, such as one parameter- 
updating attempt every one, twenty, and a hundred MC 
sweeps. They also agree with each other. Note that be- 
cause finite sizes are employed, the magnetization under 
T c at h = is also zero. With the lower parameter- 
updating frequency, the convergence was not so good and 
the error bars tend to be larger. The error bars were ob- 
tained by the jackknife method [29, 30]. These results 
suggest that the frequent parameter update does not 
make any artifacts and that it should be recommended. 

Figure 17 shows the integrated correlation time of mag- 
netization obtained at different parameter-updating fre- 
quencies. The height of data is expected to converge to 
the integrated correlation time between samples. This 
was calculated by using the jackknife method with dif- 



1.5 F 




5 -0.5 



-1 

-1.5 [ - 

1 1.5 2 2.5 3 3.5 4 4.5 5 
Temperature 

Video 1. (Color online) The red horizontal line is drawn 
between T = and 2.27 (« T c ) at h = 0. Black points 
show the sampled set of values of temperature and external 
field. Black lines are drawn when the update is accepted 
between the conditions. Parameters within 5000 MC sweeps 
are shown, during which 100 parameter-updating attempts 
were made. The windows change every 1000 MC sweeps. The 
linear lattice size L was 20. 
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FIG. 15. (Color online) Heat capacity per spin, C, at h = 0. 
The linear lattice size L was 80. As the legends shown in the 
figure, green square, blue circle, purple triangle, and light- 
blue inverse-triangle represent that one parameter-updating 
attempt is made every one, two, twenty, and a hundred MC 
sweeps, respectively. The exact result (red solid line) was 
obtained by Berg's program [30] based on Ref. [32]. 



ferent bin sizes [29, 30]. Data were stored every ten MC 
sweeps. Thus, the correlation time measured by MC 
sweep should be ten times larger. The error bars were 
obtained with the x 2 distribution. These results suggest 
that the higher parameter-updating frequency was em- 
ployed, the lower correlation time was obtained. There- 
fore, frequent parameter updates are preferred. Note that 
the observation that the frequent parameter updates are 
preferable is in accord with the statement that frequent 
replica-exchanging attempts arc recommended [33, 34]. 
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FIG. 16. (Color online) Magnetization per spin m when h = 
0. As the legends shown in the picture, the green square, 
blue circle and purple triangle represent that one parameter- 
updating attempt is made every one, twenty, and a hundred 
MC sweeps, respectively. Some error bars were slightly shifted 
horizontally to aid the eye. 
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FIG. 17. (Color online) Correlation time analysis. Error bars 
show the 95% confident interval. 



Firstly we examine the scaling behavior of the mag- 
netization. Figures 18 and 19 show the magnetization 
as functions of T and h, respectively, and we see that 
it obeys the critical behavior of m ~ \T — T C \P and 
77i ~ H 1 /* 5 , respectively. According to the scaling ap- 
proach, when Lt or L 15 / 8 /i is large enough, then the fi- 
nite effect can be negligible. In this case, Figs. 18 and 
19 imply that those conditions are given by Lt > 0.2 and 
L 15 / 8 /i > 1.1, respectively 
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FIG. 18. (Color online) Scaled m when h 
the same as used in Ref. [351. 



10 



0. The lines are 




100 1000 
il.ii 15/8 



10000 



E. Observation of crossover 

We study the crossover behavior of the phase tran- 
sitions. We calculated the magnetization by MBAR 
around the critical point. 

We employ the finite-size scaling approach, which is 
discussed in Ref. [35] . The scaling form of magnetization 
to with respect to temperature and external field is given 
by 

mL^" = ^{L l / V t, lfr + W v h) , (30) 

where t = \T — T c \/T c and L is the linear size of lattice. 
The Greek letters v and 7 stand for critical exponents. 
In the two-dimensional Ising model, (5 = 1/8, S = 15, 
v = 1, and 7 = 7/4. 



FIG. 19. (Color online) Scaled m at T = T c . 



We now study the behavior under the conditions 
slightly different from the critical point. Figure 20 shows 
the magnetization as a function of temperature near 
h = 0. As the external field was increased, the behavior 
was differentiated in the low temperature region. Even 
in the presence of weak external field, the magnetiza- 
tion obeys i 1 / 8 when the temperature was relatively high 
enough. However, with relatively strong external field, 
the scaling behavior disappears. 

Figure 21 shows the magnetization as a function of 
external field near T = T c . As the temperature was de- 
viated from T c , the behavior is differentiated in the weak 
external field region. Thus, even with slight difference 
from T c , the magnetization obeys ft, 1 / 15 when the exter- 
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FIG. 20. (Color online) Scaled m near h = 0. 



FIG. 22. (Color online) Scaled m about the critical point. 
The linear lattice size L was 160. We only display the results 
for T < T c . 
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FIG. 21. (Color online) Scaled m near T — T c 



Figure 22 illustrates the comprehensive behavior of 
(|m|) near the critical point. Note that this is a log scale 
plot. Near the h-axis (|m|) obeys [ft.] 1 / 15 and near the 
T-axis (|m|) obeys It] 1 / 8 . 

Figure 23(a) and Figure 23 (b) show the difference 
between (\m\) L 1 / 8 and 1.22(Li) 1 / 8 and that between 
(|m|) L 1 / 8 and (L 15 / 8 ^) 1 / 15 , respectively. These data 
were obtained by the 160 x 160 lattice size simulations. 
Note that the factor 1.22 comes from the exact solu- 
tion [19, 36]. According to the crossover scaling formal- 
ism [37], if t~ 15 / 8 /i is large enough, then the magneti- 
zation obeys m ~ t 1 / 8 , and if h~ 8 / 15 t is large enough 
(t~ 15 / 8 h is small enough), then it obeys m ~ ft, 1 / 15 . Fig- 
ure 23(a) shows that if the finite-size effects are negligible 
{Lt > 0.2) and t > 0.2/i 8 / 15 (i.e., t/i" 8 / 15 is large), then 
the critical behavior is m ~ t 1 / 8 . Figure 23(b) shows 
that if finite-size effects are negligible (L 15 / 8 /i 3> 0.3) 
and t < 0.2ft 8/15 (i.e., t~ 15/8 h is large), then the criti- 
cal behavior is m ~ ft 1 / 15 . Thus, Fig. 23 clearly shows 
that the line (t = 0.2/i 8 / 15 ) gives the boarder of the two 
scaling regimes. 
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FIG. 23. (Color online) Difference between magnetization 
and its expected scaling behavior about the critical point. 
The linear lattice size L was 160. (a) |mi 1/8 - 1.22( J Lt) 1/8 j is 
illustrated. The black line is t = 0.2h s/15 . The vertical gray 
line is Lt = 0.2. (b) |mi 1/8 - (L 15/S h) 1/15 \ is illustrated. 
The black line is t = 0.2h 8 ^ 15 . The horizontal gray line is 
L 15 /»/i = 0.3 



IV. CONCLUSIONS 

In this work, we introduced a two-dimensional simu- 
lated tempering in temperature and external field, which 
we refer to as Simulated Tempering and Magnetizing 
(STM). We applied it to the two-dimensional Ising model. 
During the simulations, two-dimensional random walks 
in temperature and external field were realized. The 
random walk covered a wide area of temperature and 
external field so that the STM simulations enabled us to 
study a wide area of phase diagram from a single simu- 
lation run. 

Even though the first-order phase transitions along the 
external field change did not directly occur, the tran- 
sitions happened through high temperature regions, or 
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"super critical water regions." The dimensional extension 
allowed us to overcome the difficultly of the first-order 
phase transitions. Thus, this result suggests that the di- 
mensional extension allows us to overcome the difficulty 
of crossing the first-order phase transition points with 
the ST method. The similarity between ST and REM 
implies that the dimensional extension of REM will also 
give this property (An example is shown for the case of 
a two-dimensional REM simulation in temperature and 
pressure in Ref. [3]). 

We also performed STM simulations with several dif- 
ferent parameter-updating frequencies. According to the 
convergence and sizes of error bars, the more frequent at- 
tempts should be the better choice. The calculated auto- 
correlation time also suggested that the frequent attempt 
is favorable. 

We investigated the crossover behavior of phase transi- 
tions by calculating the magnetization per spin m around 
the critical point by the reweighting techniques. The 
results showed agreement with the previous theoretical 
studies. Thus, this supports the validity of the two- 
dimensional ST method, or STM. 



(i,3) 
(i,3) 



where n and TV are the number of occupied sites and 
the total number of sites, respectively. The first term 
corresponds to the attractive energy between particles of 



With the data of the present work, we can calculate the 
two-dimensional density of states n(E,M), so that we 
can determine the weight factor for the two-dimensional 
multicanonical simulations. Therefore, we can also per- 
form the two-dimensional multicanonical simulations. 
The work is in progress. The STM method will be very 
useful for simulating spin-glass systems, and work is also 
in progress. We also remark that the present methods 
are not only useful for spin systems but also other com- 
plex systems with many degrees of freedom. It should be 
worth noting that because this method does not modify 
the energy calculation, the method should be very much 
compatible with existing package programs. 



Appendix: Lattice gas and Ising model 

The total energy of Ising model H on a square lattice 
can be converted into that of lattice gas in the following 
manner: 



(A.l) 
(A.2) 



I 

where oi = ±1 and Sj = 1, 0. If <Xj = 1, then Sj = 1 and 
vice versa. We then have 



(A.3) 
(A.4) 
(A.5) 



I 

lattice gas. The second term corresponds to the chemical 
potential of lattice gas. The last term is a constant. Here, 
we define fi = (2h — 8 J) and E g = ~4J sisj. 
Thus, free energy per spin / is given by 



I 

- ftyVj 

-l)(2 Si -l)-/ l ^(2 Si -l), 



h = -4j s * s 3 + 2J Y1 ( Si + s j) + J J2 1 ~ h Yl ( 2s * ~ 

= -4J Y S i S 3 + SJn + 2JN ~ 2hn + hN 

(i,3) 

= -4 J Y s i s j ~ ( 2h ~ 8 J)n + (h- 2J)N , 

(i,3) 
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cxp(-/3/Y) = cxp(-/3#) 

era— ±1,cti— ±l,...,fTjv— ±1 

Y exp[-f3(E g - (j,n)] exp(-/3(ft - 2J)N) 

s =1,0,si=1,0,...,sn=1,0 

= 9exp(-/3(/i-2J)iV) 

= exp(f3pN) exp(-/?(ft - 2 J) AT) , 



(A.6) 

(A.7) 

(A.8) 
(A.9) 



where p is pressure. Instead of V , N appears. The Greek 
letter stands for the Grand partition function, where 

= E ao =l,O, Sl =l,O,...,^=l,O eX PH 3 E S - HI- The last 

two equations were obtained with grand canonical en- 
sembles. Therefore, we obtain 

-f=p-(h-2J), (A.10) 
p = h- f -2J. (A.ll) 

Thus, we conclude that the canonical ensemble of Ising 
model is equivalent to the \x-T ensemble of lattice gas 
model with the following correspondence: 

p = h-f-2J, (A.12) 
/i = (2ft - 8J) , (A.13) 
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